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FIELD OF THE INVENTION 

This invention relates generally to the field of magnetic resonance imaging (MRI), 
image display, image processing and image filtering to enhance structure and reduce noise. 

BACKGROUND OF THE INVENTION 

Magnetic Resonance Imaging (MRI) has revolutionized radiological imaging of the 
intemal structures of the himian body. It has the advantage of being noninvasive and offers no 
known health hazards. A variety of MRI protocols are currently available, however image 
acquisition techniques often suffer from low signal-to-noise ratio (SNR) and/or contrast-to- 
noise ratio (CNR), Although many acquisition techniques are available to minimize these, post 
acquisition filtering is a major off-hne image processing technique commonly used to improve 
the SNR and CNR. A major drawback of filtering is that it often diffuses/blurs important 
structures along with noise. 

Two- and higher-dimensional images are currently available through sensing devices 
that operate on a wide range of frequency in the electromagnetic spectrum - from ultrasound to 
visible light to X- and y-rays (Chu et aL, Foundation of Medical Imaging , John Wiley Sons, 
Inc., New York, NY, 1993). The usefiilness of any of these modalities depends on the 
perceptability of certain features in the created images. Acquired images are often degraded by 
various types of artifacts resulting in a lowering of signal-to-noise ratio (SNR) and contrast-to- 
noise ratio (CNR). Despite high adaptivity and robustness of the human visual perceptual 
system, low SNR and CNR often degrade the information contained in the image and its utility. 
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Further, image artifacts affect many image processing tasks such as segmentation (Sahoo et al.. 
Computer Vision Graphics and Image Processing 41 :233-260 (1988); Udupa et al., Graphical 
Models and Image Processing 58(3):246-261 (1996)), registration (Wells et al. Medical Image 
Analysis 1:35-51 (1996); J. Maintz et al.. Medical Image Analysis 1:151-161 (1996)), and 
visual rendition (Udupa, Radiographics 19:783-806 (1999); Hohne et al., (eds.), 3D Imaging in 
Medicine: Algorithms. Systems, Apphcations , Springer- Verlag, Berlin, 1990), which are 
crucial in many applications. Therefore, noise reduction, that is, the improvement of SNR and 
CNR, is of considerable interest in many imaging applications. 

There are basically two approaches for improving SNR and CNR: by changing the 
method of image acquisition and by utihzing post-acquisition image processing. Improving the 
method of acquisition usually entails lengthening the acquisition time, or sacrificing spatial 
resolution, or upgrading of the imaging device hardware. Filtering is an off-line image 
processing approach, which, if properly designed, may be less expensive than, and as effective 
as, acquisition improvement approaches without having to sacrifice spatial resolution. 

The aim of filtering is to enhance wanted (structure) information and to suppress 
unwanted (noise) information. Filtering operations may be classified into two categories - 
enhancing (also know as high pass filtering), wherein wanted information is enhanced 
hopefully without affecting unwanted information, and suppressing (also known as low-pass 
fihering), wherein unwanted information is suppressed hopefixUy without affecting wanted 
information. 

Noise reducing filtering operations may be classified into two major categories. The 
first is space-invariant filtering (Rosenfeld and Kak, Digital Picture Processing , Academic 
Press, Inc., Orlando, FL, 1982), wherein a spatially independent fixed smoothing operation is 
performed over the entire image. The other is space-variant filtering (Lee, IEEE Transactions 
on Pattern Analysis and Machine Analysis 2:165-168 (1980)), wherein the smoothing operation 
is modified by local image features. 

A major drawback of space-invariant filtering techniques is that, along with noise, they 
often blur important structures. To overcome this problem, a variety of local image feature- 
dependent adaptive filtering strategies have been developed, including local image-statistics- 
driven filtering (Lee, 1980), least-squares error filtering (Chan and Lim, IEEE Transactions on 
Acoustic, Speech and Signal Processing 33:1 17-126 (1985)), gradient inverse-weighted 
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filtering (Wang et al. Computer Graphics and Image Processing 15:167-181 (1981); Wang, 
IEEE Transactions on Signal Processing 40:482-484 (1992)), multiple model Kalman filtering 
(Woods et al, IEEE Transactions on Pattern Analysis and Machine Analysis 9:245-253 
(1987)), adaptive recursive filtering using local image features (Rank et aL, IEEE Transactions 
on Image Processing 1:431-436 (1992)), local shape-based template-matched adaptive filtering 
(Ahn et aL, IEEE Transactions on Medical Imaging 18:549-556 (1999)), and gradient- 
controlled anisotropic difflisive filtering (Perona et al, IEEE Trans-actions on Pattern Analysis 
and Machine Analysis, 12:629-639 (1990)). This later method has become quite popular, and 
has been fiirther extended and utilized (Gerig et ah, IEEE Transactions on Medical Imaging 
1 1:221-232 (1992); Jackson et al, 1 Computer Assisted Tomography 17:200-205 (1993); Bajla 
et a/., Proc, MEDINFO, Part 1:683-686 (1995); Parker et al, Proc. International Society of 
Magnetic Resonance in Medicine 7:175(1999)) in a variety of applications. 

Anisotropic diffiisive fihering (Perona et aL, 1990) provides a method, in which in each 
iteration, intensity diffusion takes place between every two spatially close image spatial 
elements ('spels'). The magnitude of the diffiision between the two spels is determined by the 
magnitude of the gradient of image intensity between them. It varies in a monotonically, non- 
decreasing manner up to a certain gradient magnitude value, and subsequently changes in a 
monotonically, non-increasing fashion. The flow across boundaries is thus significantly 
reduced by utilizing information about the presence of boundaries provided by intensity 
gradients. 

A significant drawback of this very usefiil approach, however, is that it does not offer 
any image-dependent guidance for selecting the optimum gradient magnitude. More 
importantly, since it does not use any morphological or structural information to control the 
extent of diffusion in different regions, fine structures often disappear and fuzzy boundaries are 
further blurred upon filtering. 

Thus, until the present invention, there has remained a need in the art for a reliable 
method for filtering SNR and CNR in an MRI acquired image in which intensity gradients are 
accurately determined, but in which selection of the optimum gradient magnitude is permitted 
and fine structures are detected. To overcome the problems of diffusive filtering in the prior 
art, the inventors have explicitly utiUzed 'object size' information - or 'scale' - to control the 
degree of smoothing that is done in different regions of the image. The notion of scale, as 
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defined herein, is different from the term 'scale' used in scale-space computer vision literature 
(Lindeberg, Scale-Space Theory in Computer Vision , Boston, MA: Kluwer, 1994; Lifshitz et 
aL, IEEE Transactions on Pattern Analysis and Machine Intelligence 12(6):529-5407 (1990); 
McAuliffe et al, in Proc, Visualization in Biomedical Computing, Lecture Notes in Computer 
5 Science, 1131, Springer Verlag, Berlin, 1996, pp. 173-182; Pizere^ a/., Computer Visionand 
Image Understanding, 69(1):55-71 (1998)), although the intent is the same, z.e., to tailor the 
processing to local object size. 

SUMMARY OF THE INVENTION 

10 The invention provides novel space- variant anisotropic (scale-based) methods for 

filtering SNR and CNR in an MRI acquired image - a spatial-resolution adaptive scale- 
computation method, and a scale-based neighborhood averaging method, and a scale-based 
diffusive fihering method -both of the scale-based methods outperform the anisotropic 
diffiisive approach in the prior art. The present scale-based filtering methods use local 

1 5 structure size or "object scale" information to arrest smoothing around fine structures and 
across even low-gradient boundaries. One method uses a weighted average over a scale- 
dependent neighborhood, while another employs scale-dependent diffusion conductance to 
perform filtering. 

Both methods adaptively modify the degree of filtering at any image location, 
20 depending on local object scale. This permits the accurate use of a restricted homogeneity 

parameter for filtering in regions with fine details and in the vicinity of boundaries, while at the 
same time, a generous filtering parameter is used in the interiors of homogeneous regions. 

Scale-based filtering basically uses the local scale information at every spel to control 
the extent of filtering. By definition, the scale in regions with fine details or in the vicinity of 
25 boundaries is small. Thus, a restricted homogeneity parameter is used for filtering in small- 
scale regions (corresponding to fine structures and vicinity of boundaries), and to use a 
generous filtering parameter at large-scale regions (corresponding to interiors of large 
homogeneous regions). 

Consequently, the local scale information indicates the radius of the largest hyperball of 
30 a "homogeneous" region centered at the spel. This concept was originally introduced by the 
inventors at the SPIE Conference Medical Imaging 2000, and further described in Computer 
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Vision and Image Understanding, 77:145-174 (2000), where it was objectively demonstrated 
that the method significantly improve image segmentation based on fuzzy connectedness. 
Thus, the present invention provides a considerable body of evidence to show the effectiveness 
of the scale information in improving diffusive filtering, as well as for devising simple new 
5 strategies based on weighted averaging of intensities. 

Qualitative comparison experiments based on both mathematical phantoms and medical 
(patient) MR images utilizing the anisotropic diffusion method and the two scale-based 
methods, showed significant improvements using the scale-based methods over the extant 
anisotropic diffusive filtering method in preserving fine details and sharpness of object 
1 0 boundaries. Quantitative analysis on geometric phantoms generated under a range of 

conditions of blurring, noise, and background variation confirmed the superiority of scale- 
based approaches embodied in the invention. 

Thus, a scale-based filtering method of the invention is provided, wherein an enhanced 
MRI-acquired image is achieved for any selected MRI protocol, and for any selected region of 
a patient's body, by the steps comprising: imaging the region of the patient's body by the 
selected MR protocol to form an acquired image; and filtering the acquired image by a scale- 
based spatial resolution adaptive method using region-constancy based on local homogeneity to 
produce an enhanced image. The resulting image is enhanced in a manner independent of 
variations within or between patients, within or between tissues being imaged, or within or 
between MR devices used to acquire the image. 

Also provided are computerized systems and devices for implementing the foregoing 
methods. 

In addition, as provided in at least one embodiment, scale-based filtering methods of the 
invention are automated. Moreover, an embodiment of the invention is provided in which 
scale-based filtering methods are built into an MR scanner, permitting production of enhanced 
real time images. 

Also provided is an enhanced MR image for imaging a region of the body of a patient 
according to one or more of the embodied scale-based filtering methods of the invention, 
including a scale-based neighborhood averaging method, or a scale-based diffusive filtering 
method. 
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Further provided are apparatus or devices incorporating one or more of the foregoing 
features, including one or more of the embodied scale-based filtering methods of the invention 
to produce an MR image in which fine details and sharpness of object boundaries are preserved 
under a variety of conditions. 
5 Additional objects, advantages and novel features of the invention will be set forth in 

part in the description, examples and figures which follow, and in part will become apparent to 
those skilled in the art on examination of the following, or may be learned by practice of the 
invention. 

1 0 DESCRIPTION OF THE DRAWINGS 

The foregoing summary, as well as the following detailed description of the invention, 
will be better understood when read in conjunction with the appended drawings. For the 
purpose of illustrating the invention, there are shown in the drawings, certain embodiment(s) 
which are presently preferred. It should be understood, however, that the invention is not 

15 limited to the precise arrangements and instrumentaUties shown. 

FIG. 1(a) and 1(b) depict local structure size or scale. FIG. 1(a) depicts a 2D slice from 
a MR 3D scene of the head of a multiple sclerosis patient. Using a suitable region- 
homogeneity parameter, a local determination is done to determine the largest disc that can be 
centered at any point, such as p, within which the intensities are homogeneous. The scale at p 

20 is bigger than that at q or at o. Scale is usually small in the vicinity of boundaries, such as at o, 
FIG. 1(b) depicts the scale scene of the MR slice shown in FIG. 1(a). In FIG. 1(b), intensity at 
a location is proportional to the scale value at that location. 

FIG. 2(a)-2(f) depict diffusion flow across blurred object boundaries and the usefiihiess 
of object scale information to arrest such flows. FIG. 2(a) depicts an initial binary phantom 

25 scene. FIG. 2(b) depicts a test phantom scene after blurring and adding a zero mean Gaussian 
noise and a slow varying (ramp) background component across columns. FIGs. 2(c) and 2(d) 
depict the scenes that resulted firom filtering the scene in FIG. 2(b) by anisotropic diffusion for 
3 and 10 iterations, respectively. FIGs. 2(e) and 2(f) depict the scenes that resulted from 
filtering the scene in FIG. 2(b) by scale-based diffiision for 3 and 10 iterations, respectively. 

30 FIG. 3(a) graphically depicts the shape of diffusion conductance function that complies 

with the principle of nonlinear diffusion, FIG. 3(b) graphically depicts the shape of diffusion 
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flow magnitude function that complies with the principle of nonlinear diffusion. FIG. 3(c) 
graphically depicts the shape of different conductance functions at different scales. The right- 
most curve of FIG. 3(c) represents the conductance function at maximum scale tmax- FIG. 3(d) 
graphically depicts the shape of different diffusion flow magnitude functions at different scales. 
5 The right-most curve of FIG. 3(d) represents the diffusion flow magnitude at maximum scale 

^MAX. 

FIGs. 4(a)-4(h) depict comparisons among the three filtering methods using the 2D slice 
of a MR 3D scene of the head of a multiple sclerosis patient originally shown in FIG. 1(a). 
FIGs. 4(a) and 4(b) depict 2D scenes within two ROIs, one from the top left comer and the 

10 other from the bottom right comer of the original scene in FIG. 1(a), FIGs. 4(c)-4(h) depict 
filtered scenes (in 2D) corresponding to the scenes in FIGs. 4(a) and 4(b) that resulted from 
anisotropic diffusion (shovra in FIGs. 4(c) and 4(d)), from scale-based neighborhood averaging 
(shown in FIGs. 4(e) and 4(f)), and from scale-based diffusion (shown in FIGs. 4(g) and 4(h)), 
FIGs. 5(a)-5(e) depict comparisons among the three filtering methods using a 3D MR 

1 5 scene of the head of a multiple sclerosis patient, FIG. 5(a) depicts a 2D slice of the original 3D 
MR scene and the selected ROI, FIG. 5(b) depicts the scene within the ROI, magnified. FIGs. 
5(c)-5(e) depict filtered scenes corresponding to the scene in FIG, 5(b) that resulted from 3D 
filtering using anisotropic diffusion (shown in FIG, 5(c)), using scale-based neighborhood 
averaging (shown in FIG. 5(d)), and using scale-based diffusion (shown in FIG. 5(e)). 

20 FIGs. 6(a)-6(f) depict comparisons among the three filtering methods using a 2D 

phantom scene containing different geometric shapes of different size and orientation, FIG. 
6(a) depicts an original phantom scene (one among 25 such scenes) at a low level of blurring 
and a medium level of noise, together with an intensity profile along the depicted bright 
horizontal line. Six ROIs are also shown within FIG, 6(a), which are used for FIG. 7. FIG. 

25 6(b) depicts a magnified version of a section of the profile in FIG, 6(a), which section is shown 
at the bottom of FIG. 6(a) as a rectangle outhned in black. FIG. 6(c) depicts the profile 
corresponding to that in FIG. 6(b), but which were obtained for the original binary phantom 
scene after adding only the background ramp component (but not adding blurring and noise). 
FIG. 6(d)-6(f) depict profiles corresponding to that in FIG. 6(b), but which were obtained from 

30 the filtered scenes: by anisotropic diffusion (shovra in FIG. 6(d)), by scale-based neighborhood 
averaging (shown in FIG. 6(e)), and by scale-based diffusion (shown in FIG. 6(f)). 
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FIGs. 7(a)-7(d) depict comparisons among the three filtering methods using the 2D 
phantom scene of FIG. 6(a) for the six ROIs shown therein. FIG. 7(a) depicts scenes within the 
six ROIs taken fi-om the unfiltered phantom scene. FIGs. 7(b)-7(d) depict the corresponding 
filtered scenes that resulted fi:om anisotropic diffusion ((shown in FIG. 7(b)), fi-om scale-based 
neighborhood averaging (shown in FIG. 7(c)), and fi*om scale-based diffiision (shown in FIG. 
7(d)). 



DESCRIPTION OF PREFERRED EMBODIMENTS OF THE INVENTION 

The invention provides novel scale-based MR image filtering methods, including scale- 
based neighborhood averaging and scale-based diffiision, which avoid or prevent blurring 
around fine structures and across boundaries using local structure size or "scale" information. 
Fundamentally, local structure size, or "scale", as defined by region homogeneity, improves the 
filtering operation by its utility in controlling the operation. The scale concept allows a 
restricted smoothing selectively to be achieved in the vicinity of boundaries and in fine 
structured regions, while at the same time allowing more generous filtering in the interiors. 
Thus, a scale-based filtering method of the invention is provided, wherein an enhanced MRI- 
acquired image is achieved for any selected MRI protocol, and for any selected region of a 
patient's body, by the steps comprising: imaging the region of the patient's body by the 
selected MR protocol to fomi an acquired image; and filtering the acquired image by a scale- 
based spatial resolution adaptive method using region-constancy based on local homogeneity to 
produce an enhanced image. 

In the discussion which follows, both scale computation and the filtering operation are 
set up in a general n-dimensional firamework, which also takes into account the anisotropy of 
spatial resolution. The fi-amework is designed in such a way that there is only one parameter to 
be specified - the sole image-dependent homogeneity parameter - and methods are provided 
for selecting this parameter utiUzing either user guidance or global image statistics. However, 
to better understand the present disclosure, certain terms and notations are defined as follows. 

Most digital imaging systems produce images on a rectangular grid where the basic 
image elements, referred to as "spels," are hypercuboids. The term "spels" is an abbreviation 
for spatial elements; in 2D spels are 'pixels,' while in 3D they are 'voxels.' Mathematically, 
this is explained by letting w-dimensional EucUdean space 9?" be divided into hyper cuboids by 
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n orthogonal families, each of equally distant parallel planes. Let vbe an «-tuple whose 
element Vi represents the distance between two successive parallel planes perpendicular to the 
coordinate axis. (A 'tuple' is a recognized term, simply referring to a mathematical relation 
between its attributes.) In other words, n is the size of each spel along the direction of the 
coordinate axis. 

The term vas used herein means the resolution of the grid system. Further, a coordinate 
system is chosen such that the centroid of each spel has the coordinates of the form (ci vi, C2V2, . 
. .,Cn v„) where c,'s are integers. Thus, a one-to-one and onto mapping between the set of spels 
and as follows: an /j-tuple c, denoting a point in 2", is mapped to a spel, such that cj Vj, for 1 
<j < n, gives they* coordinate of the center of the spel represented by c. With this 
interpretation of spels, Z" itself represents the set of all spels of 9?", and concepts of spels and 
points denoting their centroids may, therefore, be used interchangeably to describe the 
invention. 

For computational simplicity, a hard relation of adjacency was used for a in all 
experiments; however, fuzzy relations for a should also be considered, especially those 
functional forms for which closely resemble the point-spread function of the imagmg 
device. A fuzzy relation a on is said to be a "fuzzy spel adjacency" if it is reflexive and 
symmetric (for a description of 'fuzzy relations,' see Kaufmann, Introduction to the Theory of 
Fuzzy Subsets , 1, Academic Press, New York, 1975). It is desirable that a be such that its 
membership value judc, d), for dny c, d & Z", is a non-increasing function of the distance \\c-d\\ 
between c and d, where || • || denotes any L2 norm in 91". 

The functional form of //q, that is used for the image processing results presented herein, 
is as follows (however, a general functional form for is considered for all formulations). For 
any spel c, /udc, c) = 1 . Further, for any spels c and d, ndc d) = \, '\ic2sAd differ in exactly 
one coordinate by 1. Otherwise, //^(c, tT) = 0. 



In 2D: 




1, if c = J or c, J are 4-adjacent, 
0, otherwise. 



(Equation 1) 
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Thus, it is readily seen that a is indeed a fuzzy spel adjacency. The triplet (Z^, a, v) is, 
therefore, called a "fuzzy digital space," in which a is any fuzzy spel adjacency, and vis the 
resolution of the grid system. A fuzzy digital space is a concept that characterizes the 
underlying digital grid system, independently of any image-related notions. 

5 A "scene over a fuzzy digital space" (Z", a, v) is a pair C = (C,/), where C= {c | -bj 

< cj < bj for some g Z" }, Z" is the set of w-tuples of positive integers, /is a function whose 
domain is C, called the "scene domain", and whose range is a set of integers [i, //]. The terms 
"scene" and "image" are used interchangeably, herein, as recognized in the art. The invention 
generalizes in a straightforward manner to vector-valued scenes, z.e., when/is a vector- valued 

1 0 function. However, for the sake of simplicity of description,/is treated herein as a scalar- 
valued function. A scene ^over the triple (Z", a, v) is referred to as a "binary scene" if the 
range of /is {0, 1}. 
1) Anisotropic Diffusive Filtering. 

A diffusive filtering process can be defined using the divergence operator 'div' on a 

15 vector field, A mathematical formulation of the diffusion process on a vector field V at a point 
c in coordinate-fi-ee form is given by: 

^ = div V = lim f V • ds, (Equation 3) 

20 wherein Ar is the volume enclosed by the surface s (surrounding c\ and d% = \xds, wherein u 
is the unit-vector normal and outward-directed to the 'infinitesimal' surface element ds. The 
diffusion process is controlled by the intensity flow vector field V. The functional form for 
flow is defined as 

25 V = GF, (Equation 4) 

wherein G is the diffusion-conductance function, and F is the scene intensity gradient vector 
field. 

In an isotropic diffusion, the conductance function G is constant, and it has been argued 
30 by Perona et al, 1990, that such diffusion methods blur object boundaries and structures. As a 
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result, Perona et al, proposed an anisotropic diffusion method in which G is spatially varied as 
a nonlinear function of the magnitude of the scene intensity gradient, such that the conductance 
function encourages smoothing within a region (characterized by low intensity gradients) in 
preference to smoothing across boundaries (characterized by high intensity gradients). 
5 2) Scale-Based Filtering, 

Although, the main motivation for anisotropic diffusive filtering was to reduce noise 
while minimizing blurring across boundaries, since it does not account for local structure size, 
it often loses/blurs fine structures and diffuses across boundaries. This phenomenon is 
illustrated in FIGs. 2, 4-7, which will be discussed in greater detail below. 

10 In this section, however, two scale-based filtering methods are presented that use the 

structure size information to accurately arrest diffusion around fine structures, and across even 
low gradient boundaries. The size of local structures is determined at every spel under a pre- 
specified scene-dependent, region-homogeneity criterion, and then the generosity of filtering is 
controlled using this parameter. For example, m FIG. 1(a), which is a 2D slice from a MR 3D 

15 scene of the head of a multiple sclerosis patient, using a suitable region-ho|mogeneity 

parameter, a local determination is done to determine the largest disc that can be centered at 
any point, such as /?, within which the intensities are homogeneous. In FIG. 1(a), the size of the 
local structure to which pixel p belongs is bigger than that to which qoxo belong. The scale is 
usually small in the vicinity of boundaries, such as at o. It is, therefore, reasonable to allow a 

20 more generous filtering at p (i.e., an averaging over a larger neighborhood, or a diffusion using 
a higher conductance), than that at q or at 

The theory of defining scale as the size of the local structure and utilizing scale in 
carrying out local operations was originally proposed by the inventors (Saha et al, 2000, herein 
incorporated by reference). Following the same line of reasoning, "scale" in a scene ^at any 

25 spel, CG C, is defined as the radius, r(c), of the largest hyperball centered at c, which lies 

entirely in the same object region (defined under a pre-specified region-homogeneity criterion). 
Ironically, it appears that object definition is first needed to determine scale. Thus, a simple 
and effective algorithm (presented in the following section, as Algorithm OSE) has been 
defined by the inventors (see also Saha et al, 2000), which estimates r{c) at every spel, cgC, in 

30 any scene C without expUcit object segmentation, but based on continuity of region 

homogeneity. This 'rough' object size information is utilized to control the diffusion process 
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SO as to improve diffusive filtering and to devise a simple neighborhood averaging strategy that 
can also out-perform diffusive filtering. 
2.1) Computation of Scale. 

A "hyper ball Bk, Jic) of radius A: > 0 and with center at c e C in a scene C=(Qf} over 
(Z^ a, v) is defined by: 



eeC 



(Equation 5) 



wherein the term e is defined within the equation as a spel belonging to C. 

Note that B],^Xc) forms a digital hyper ball when the spatial resolution is isotropic, and a 
10 digital hyper ellipsoid when it is not. The equation is formulated so that, in both cases, 5k,v(c) 
represents a (quantized) hyper ball in the scene domain. For the purposes of the invention, 
Bk,^ic) is referred to as a hyper ball, irrespective of whether or not vis isotropic. 

For a hyper ball Bk^J^c) of any radius k>0, and centered at c, a fraction, FOk^^c) 
("fraction of object"), is defined to indicate the fraction of the ball boundary occupied by a 
1 5 region which is sufficiently homogeneous with c, by: 



FOf^^ (c) = — ^ — — j ^ (Equation 6) 



5,/c)-5,_,,(c)| 



wherein 1 5^, ^c) - Bk.\,Jic) \ is the number of spels in {i.e., the volume of) Bk, Jl^c) - Bk-iA^X 
20 and Wi^ is a homogeneity function. A detailed description of the characteristics of is found 
in Saha et al, 2000, herein incorporated by reference. The membership functions utilized in 
fuzzy subset theory corresponding to the fuzzy proposition "x is small" can be used for W^. A 
zero-mean, unnormalized, Gaussian function is a preferred choice for Wy^. 

The algorithm for object scale estimation is presented below as Algorithm OSE, 
25 wherein "OSE" refers to "Object Scale Estimation." It has been modified firom that which was 
originally presented by Saha et al, 2000, to account for resolution anisotropy of acquired 
scenes. 

Algorithm OSE 
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Input: C,ceQWif,,di fixed threshold ts. 

Output: r(c). 

begin 

SQtk= 1; 

5 while FOk, v(c) > ts do 

set A:to A:+ 1; 
endwhile\ 
set r{c) tok-\\ 
output 
10 end. 

The algorithm iteratively increases the ball radius A: by 1, starting from 1, and checks for 
FOk, v(c), the fraction of the object containing c that is contained in the ball. The first time that 
this fraction falls below ts, the ball then contains an object region different from that to which c 
belongs. 

15 In accordance with Saha et al, 2000, ts = 0.85 is used. The rationale for this choice is 

as follows: In a 3x3 neighborhood of a pixel c, in a 2D scene, one out of the eight 
neighboring pixels may belong to a different object (to account for noise), but the 
neighborhood is still considered to be entirely within the same object. This leads to ts = %. No 
significant change in the results was observed by varying ts in the range 0,8 ^ ^ 0.9. 

20 To complete the description of scale computation, the following is a description of two 

methods by which a proper value is selected for region-homogeneity parameter cr^. In the first, 
which is an operator interactive training-based method, an operator using a mouse controlled 
brush, paints on a display of a slice of an image to specify objects of interest, excluding inter- 
object boundaries. The mean, Mh, and the standard deviation, ah, of local intensity gradient 

25 magnitudes, are computed as the mean and the standard deviation of the intensity differences, 
\A^) -Ad) I , for all possible pairs (c, d) of spels in the painted regions, such that c^dmd 
jLia{c,d) > 0. These estimates are then used to set up the value of the homogeneity parameter as 
follows: 

30 ay, = Mh + thCTh, (Equation 7) 
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wherein 4 is a constant, and its value is set to 3. This is based on the rationale that in a normal 
distribution, 3 standard deviations on the two sides of the mean cover 99.7% of the population. 

The second method of selecting cr^ is similar to the preceding method, except that it 
requires no operator intervention. Instead of computing the mean and the standard deviation of 
local intensity gradients over the painted regions, the local, in-plane intensity gradients (z.e., c 
and d are in the same plane) are first determined over the entire scene domain. Local, out-of- 
plane intensity gradients are excluded because, in most medical imaging applications, the out- 
of-plane resolution is relatively coarse. However, local out-of-plane intensity gradients may be 
used when the spatial resolution is isotropic. The upper 10th percentile of the gradients are 
then discarded in order to account for inter-object boundaries. The mean and the standard 
deviation are computed over the lower 90 percentile of the gradients. Finally, the expression in 
Equation 7 is used for determining <jy,. In this fashion, the only parameter required in the 
method is determined automatically. 

In the examples that follow, the second method was used for computing a^. For 
example, the scale image shown in FIG. 1(b) corresponds to the MRI sUce of the head of the 
patient shown in FIG. 1(a). In the scene in FIG. 1(b), intensity at a location is proportional to 
the scale value at that location. As can be seen, the interior regions have higher scale (brighter 
intensities), while in contrast, the detailed structural regions, as well as the regions in the 
vicinity of boundaries, have lower scale (darker intensities). 

2.2) Scale-Based Neighborhood Averaging. 

A fixed (Gaussian) kernel-based filtering method essentially convolves the original 
scene with a fixed kemel. A major drawback of such a method is that, when the kernel is 
defined over a small neighborhood, the method fails to significantly increase the SNR upon 
filtering. On the other hand, when the kemel is defined over a large neighborhood, smaller 
shapes are lost and substantial amount of blurring takes place across boundaries. Finding an 
optimum neighborhood size (kemel) had been a fimdamental problem in this type of approach 
in the prior art. By bringing object scale into the fi-amework, not only is the size of the 
neighborhood specified, but also it is allowed to vary over the scene domain in a meaningful 
way. This allows large neighborhoods to be used in the interiors of large structures (e.^., at/? in 
FIG. 1(a)), while small neighborhoods are used in fine, detailed regions and in the vicinity of 
boundaries {e.g., at q or at o) in the same figure. 
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A scale-based neighborhood averaging method can be formally defined as creating an 
output filtered scene C = (C,/sa.) for an input scene {C,f) over (Z«, a,v\ where, for any 

C 6 C, 



wherein coc is a distance-dependent weighting fimction. The true distance between c and e is 
used for the weighting function cOc, with the purpose being to develop a scale-dependent 
filtering kernel using 5^c),v{c). The filtering kernel cOc should satisfy the following properties: 

1 . The domain of 6)^, is the set of non-negative numbers. 

2. The range of cOc, is [0,1] and codO) = 1. 

3. There is one control parameter of q)c, such that for a given distance between c and 
e, COc, is monotonically non decreasing with . 

4. For any given control parameter a^^ , cOc is monotonically non-increasing with the 
distance between c and e. 

In the fuzzy subset theory, the membership functions corresponding to the fuzzy 
proposition "x is small" satisfy these properties. Three such popular functions are (1) 
unnormaUzed Gaussian, (2) ramp, and (3) box window. In the examples that follow, an 
unnormaUzed Gaussian function for cOc was used. Intuitively, it follows that should be a 
function of r(e), and accordingly cr^^ = r{c) is used. 

2.3) Scale-Based Diffusive Filtering. 

In the original anisotropic diffusion method, flow across objects is reduced due to high 
intensity gradients at boundaries. Intensity gradients at spels on blurred boundaries, where 
gradient is distributed over a thick boundary band, are often not high enough to arrest flow. As 
a result, diffusion flow takes place across blurred boundaries. The effect of this phenomenon 
becomes visually prominent when the filtering method is executed for several iterations. 
Witiiin an 'iteration,' estimated intensity flow diffuses between every two adjacent spel in the 




(Equation 8) 
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image domain. This effect is demonstrated in FIG. 2(d), wherein the diffusion process is run 
for 10 iterations on the 2D scene shown in Figure 2(b). 

Thin and fine structures in images often fade or disappear, even after fewer than 10 
iterations of anisotropic diffusion. The behavior of anisotropic diffusion around thin and fine 
structured regions is further discussed in the Examples that follow; see also FIGs. 4, 5 and 7. 

The purpose behind scale-based diffusion is to effect a restricted diffusion around thin 
and fine structured regions and in the vicinity of boundaries, while at the same time allowing a 
generous diffusion in the interiors of homogeneous regions. This is achieved in the present 
embodiments by using a variable diffusion conductance function that is controlled by the scale 
parameter. As demonstrated in the section entitled Scale Based Filtering and shown in FIG. 1, 
scales in the vicinity of boundaries and in thin and fine structured regions are small, while with 
the interiors they are large, thus showing the effectiveness of a scale-dependent diffusion 
process. Compare FIG. 2(c) with 2(e), which were obtained after 3 iterations, and FIG. 2(d) 
with 2(f), which were obtained after 10 iterations. See also FIGs. 4, 5, and 7. 

However, before further describing the scale-based method, a discussion of the 
characteristics of diffusion conductance flmctions is needed. 

Perona et al, 1990, described qualitatively the shapes of the diffusion conductance 
function, G, and diffusion flow magnitude fimction, |V| (FIGs. 4 and 6 in Perona et al, 1990) 
that comply with the idea behind nonlinear diffusion. The required shapes of diffusion 
conductance and flow magnitude functions are shown in FIGs. 3(a) and (b), respectively. The 
necessary properties of G may be formulated as follows: 

1 . The domain of G is the set of all non-negative numbers. 

2. The range of G is [0, 1], and G — 1 when |F| = 0, while G 0 when |F| — > oo. 

3. There is a control parameter cr, such that for any given |F| and o\ < cPi, G^ <G^ , 



Further, |V| (see Equation 4) has a maximum value at |F| = cr, and |V| is 
monotonically increasing for |F| < a; while it is monotonically decreasing for |F| >cr. 



Intuitively, it appears that G should be monotonically decreasing with |F|, However, not 
every monotonically decreasing function results in a valid flow function. For example, if G is 



selected in the expression G = 




the resulting conductance function is monotonically 
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decreasing, but it does not result in a valid flow function because the latter is a throughout- 
increasing function of |F|, and it possesses no finite optimum. Therefore, care must be taken in 
selecting a conductance function so that the resulting flow function possesses a unique 
maximum. Two valid conductance functions were proposed in Perona et al, 1990, and Gerig et 
5 al, 1992. 



10 



-— — TTg- |^>0, (Equation 9) 



1 + 



Fr 



and G2 (F) = e . (Equation 1 0) 



Both functional forms in both anisotropic and scale-based diffusion methods have been 
applied to patient MR images, and no significant visual differences were produced in the 
results. Consequently, the unnormahzed Gaussian function of Equation 10 was used herein for 
all evaluations. 

1 5 As described above, cris the control parameter that determines the generosity of 

filtering. When cris large, the generosity of filtering is high, and the possibihties of blurring 
across boundaries increase. On the other hand, when cris small, the generosity of filtering is 
low, and more noise survives after filtering. Following this principle, crmay also be thought of 
as a region-homogeneity parameter. 

20 As previously pointed out, the basic premise behind scale-based diffusion is to allow 

only a restricted diffusion at small-scale regions, while at the same time being generous at high- 
scale regions. This is achieved by using an adaptive homogeneity parameter oj, in place of crin 
Equation 10, as a monotonically, non-decreasing function of local scale. Consequenfly, was 
used herein as a linear function of scale. 

25 Given this understanding of scale-based diffusion conductance, the process of scale- 

based diffusion is completely described by Equations 3, 4, and 10. Accordingly, throughout the 
remainder of the disclosure, the details of the scale-based diffusion process are presented as 
applicable to a scene over (Z„, a, v). Like anisotropic diffusion, scale-based diffusion is also an 
iterative process, and let t denote the iteration number. Let (C,f) be the given scene, and 
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5 



let Q = {C,f^ denote the scene resulting from the diffusion process at the t^^ iteration. For any 
c,d e C, such that c^dmd jUq(c, d)^0, let D(c, d) denote the unit vector along the direction 
from the point c toward the point d. Let F/(c, d) represent the intensity gradient vector along 
D(c, d) given by 

^Xc,d) = -^IMzIi^=\i{c,d) (Equation 11) 

^^^•-min,[vf] 

Note that the preceding expression takes into account resolution anisotropy. 

While determining the intensity flow from spel d to spel c, the effective scale, reff{c, d), 
10 defined as follows, was utilized. 

reff{c, d) = min [r(c), r(d), r(e)] (Equation 12) 

wherein e e C is the neighboring spel of c, just opposite to d, defined hye = c-(d-c). When 
15 c is in the border of the scene domain, e was let equal to c. Using this definition of effective 
scale, the diffusion conductance function, Gt{c, d) for the flow from d to cat the iteration, is 
defined by: 



20 



25 



G,(c,J)-e'^^^^^''^> (Equation 13) 

wherein oi(c, d) is the scale-adaptive region-homogeneity parameter used for controlling the 
intensity flow from d to c, given by 



a^c.d) = ^ (Equation 14) 



In this case ruAx is an upper limit on the scale in C, and in the following examples rMAX = 5. 
cr^ is the homogeneity parameter defined above in Section 2.1, entitled Computation of Scale . 
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Accordingly, intensity flow vector (c, d) from spel c to spel d at the iteration is 
defined by: 

V,(c, d) = F,(c, d) Gt (c, d) (Equation 15) 

5 Note that (c, J)-D(c, d) is positive when ft(c)>ft (d); it is negative otherwise; and it is zero 
when c = dor ft(c) =ft(d). Thus, the flow direction between two spels is always such that it 
tries to reduce the gradient between them. FIG. 3(c) shows diffusion conductance functions at 
different scales while FIG. 3(d) shows the magnitude of intensity flow vector at different 
scales. As shown in the two figures, at the maximum scale ruAXy conductance is highest and the 
10 flow magnitude is maximum. Moreover, both diffusion conductance and flow magnitude 
decrease as scale decreases. 

With these formulations, the scale-based iterative process is defined by: 



15 



20 



f/(c), for^^O, 



where Kd is a constant. It can be shown in a manner analogous to the proof in Gerig et al, 
1992, that in order to guarantee a "monotonic variation of spel intensities" (/.e., the intensity of 
any spel ceC is either non-increasing or non-decreasing with iteration t, but is not oscillatory), 
the 'upper limit' oiKo should satisfy the following inequality. For any ce Z", 



Kj^ < = ^— — - . (Equation 1 7) 

Although Kd has no lower hmit, more iterations are needed to achieve the same degree 
of filtering when a small value for is used for Kd^ In the following examples, the value Kb was 
25 used at its upper limit. Thus, according to the adjacency criterion (see Equations 1 and 2), Kd = 
V5 in 2D m&KD^^hm 3D. 

The embodied methods of the invention are useful when implemented by recognized 
techniques to produce an automated, scale-based filtered image, and are particularly useful 
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when installed by recognized methods in the art into an MR scanner, permitting production of 
enhanced real time images. 

Also provided in the invention is an enhanced MR image for imaging a region of the 
body of a patient accordmg to one or more of the embodied scale-based filtering methods of the 
invention, including a scale-based neighborhood averaging method, or a scale-based diffusive 
filtering method. Further provided in the invention are computerized apparatus and devices 
incorporating one or more of the embodiments of the invention, including one or more of the 
embodied scale-based filtering methods of the invention to produce an MR image in which fine 
details and sharpness of object boundaries are preserved under a variety of conditions. 

EXAMPLES 

The present invention is further described in the following examples in which both 
qualitative and quantitative experiments were conducted to assess the relative performance 
5 among the three methods - (1) anisotropic diffusion, (2) scale-based neighborhood averaging, 
and (3) scale-based diffusion. These examples are provided for purposes of illustration to those 
skilled in the art, and are not intended to be limiting unless otherwise specified. Moreover, 
these examples are not to be construed as limiting the scope of the appended claims. Thus, the 
invention should in no way be construed as being limited to the following examples, but rather, 

10 should be construed to encompass any and all variations which become evident as a result of 
the teaching provided herein. 

The principles and algorithms for both scale-based neighborhood averaging and scale- 
based diffusion were as described above. Also as described, were how to determine the sole 
image-dependent parameter, Oy,, and how to fix other image-independent terms. For the 

15 anisotropic diffusion method, essentially the method of Gerig et al, 1992 was followed, except 
that in computing the image gradients, resolution anisotropy was considered as in Equation 11. 

All experimental results presented in this paper were obtained using the functional form 
for diffusion conductance of scale-based neighborhood averaging in Gerig et al, 1992, which 
was the same as Equation 10, above. Since the Gerig et al method did not describe selection 

20 of parameter at, the following was used: k= ^ll ct^^ 42 (jy^m Equation 10 - the same 

parameter that was used for scale-based diffusion. Except for the result in FIG. 2(d), all results 
for anisotropic diffusion were obtained by applying the method for 3 iterations as 
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recommended in Gerig et aL, 1992, which is herein incorporated by reference. Except for the 
result in FIG. 2{f), all results for scale-based diffusion were obtained by applying the method 
for 7 iterations. 

In Example 1, the three methods are qualitatively compared using a phantom image and 
5 MR images of the heads of multiple sclerosis patients. A quantitative comparison based on a 
phantom containing different geometric shapes is presented in Example 2, to more objectively 
assess the relative performance of each of the three methods. 

Example 1 - Quahtative Comparisons. 
10 In this example, two case studies are presented from MR patient studies, and then two 

phantom experiments are described. 
Comparative Case Study 1.1. 

In FIG. 4, the filtered output from the three methods is illustrated for two regions of 
interest (ROI) - one from the top-left comer and the other from the bottom-right comer - 

1 5 selected from the MR slice displayed in FIG. 1(a). The two ROIs, and the results from the 

three filtering methods, are shown in FIGs. 4(a)-4(h). FIGs. 4(c)-4(h) depict filtered scenes (in 
2D) corresponding to the scenes in FIGs. 4(a) and 4(b) resulting from anisotropic diffusion 
(shown in FIGs. 4(c) and 4(d)), from scale-based neighborhood averaging (shown in FIGs. 4(e) 
and 4(f)), and from scale-based diffusion (shown in FIGs. 4(g) and 4(h)). As can be seen from 

20 all six filtered images (FIGs. 4(c)-4(h)), the smoothness in homogeneous regions seems 

visually to be the same in all images, while boundary contrast seems to be higher using scale- 
based methods, as opposed to using anisotropic diffusion. Also, scale-based methods were 
seen to preserve fine details (for example, the linear bright structure in the left, middle part of 
FIG. 4(b)) better than anisotropic diffusion. Moreover, among the two scale-based methods, 

25 differences did not appear to be visually significant. 
Comparitive Case Study 1.2. 

FIG. 5 illustrates SD-comparative results among the three methods. All methods were 
applied on a proton density-weighted 3D MR scene of the head of a multiple sclerosis patient. 
The out-of-plane resolution was seen to be almost 3.5 times more coarse than the in-plane 
30 resolution. The results, thus, also demonstrate the behavior of the three methods under 
resolution anisotropy. 
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The size of the scene domain was 181x246x55, with a voxel size of 0.86x0.86x3.0 
mm^ In the non-optimized experimental implementation, the computation of scale took 10-12 
seconds on a 450 MHz Pentium PC for the scene in FIG. 5, while each iteration of anisotropic 
and scale-based diffusion required 7-8 seconds. Compared to the anisotropic diffusion method, 
5 scale computation was clearly an extra step for the scale-based diffusion method. For the same 
scene, scale-based neighborhood averaging method took 17-18 seconds in total on the same 
machine, which is faster than anisotropic diffusion, which required 21-24 seconds when run for 
three (3) iterations. Moreover, the same improvement of SNR values was achieved in 
homogeneous regions whether anisotropic or scale-based diffusion was applied for the same 

10 number of iterations. However, more blurring was seen across boundaries and in fine 
structured regions in the anisotropic diffusion method. 

Since it was not clear how to best display such 3D scenes (short of showing all slices at 
sufficiently high magnification), one slice is shown in FIG. 5(a) taken from the original scene 
with an ROI. A portion of the sUce was zoomed up (magnified) in FIG. 5(b), to allow a closer 

15 scrutiny of the images. FIGs. 5(c)-(e) show the same ROI for the filtered scenes obtained using 
anisotropic diffusion (FIG. 5(c)), using scale-based neighborhood averaging (FIG. 5(d)), and 
using scale-based diffusion (FIG. 5(e)) methods, respectively. While there are no visually 
significant differences among the three methods in the level of smoothing within homogeneous 
regions, it can been seen that boimdary contrast is significantly improved using the scale-based 

20 methods. Some of the lesions (hyper intense blobs) appear blurred in the filtered scene as a 

result of anisotropic diffusion, while by comparison, their contrast is retained using scale-based 
methods. 

Although there is major resolution by anisotropy, the scale-based methods were able to 
retain boundary contrast and subtle structures. As between the two scale-based methods, the 
25 result using scale-based diffusion appears visually more pleasing than the result produced using 
scale-based neighborhood averaging. 

Phantom Example 1 . 1 . 

FIG. 2 demonstrates how diffusion proceeds across certain blurred object boundaries in 
the method of anisotropic diffusion (FIGs. 2(c) and 2(d)), whereas it is arrested in the scale- 
30 based method (FIGs. 2(e) and 2(f)). The 2D scene, shown in FIG. 2(b) was a test geometric 
phantom scene created after blurring and adding a zero mean Gaussian noise and a slow 
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varying (ramp) background component across columns to imitate white matter and gray matter 
regions in the human brain. One sUce was selected from a 3D proton density MR scene of the 
head of a multiple sclerosis patient undergoing routine clinical evaluation. Subsequently, the 
white matter and gray matter regions in the original scene were carefully segmented using a 
5 user-steered segmentation method (Falcao et aL, Graphical Models Image Processing 60:233- 
260 (1998)). To the white matter region in the scene, a value was assigned equal to the average 
of the intensities within the white matter region in the original MR slice, and to the background 
region within the brain, a value was assigned equal to the average of the intensities within the 
gray matter region in the original MR shce. 

10 From this starting scene, the 2D phantom scene was generated by applying a Gaussian 

blurring and a zero-mean Gaussian noise, and finally by adding a fixed slow varying (ramp) 
background component across columns. FIGs. 2(a) and 2(b) represent the 2D scenes before 
and after adding noise, blurring, and background variation, respectively. The size of the scene 
was 185x232 with an in-plane resolution of 0.86x0.86 mm^. FIGs. 2(d) and 2(f) were obtained 

15 by applying the anisotropic diffiision and scale-based diffusion methods, respectively. 

In both cases, 10 iterations were used. A few more iterations than the recommended 
number were used for the anisotropic method to emphasize the diffusion across boundaries. On 
the other hand, when the process was run for fewer iterations, /.e., the recommended number of 
3 iterations, the noise reduction was not significant, as can be seen from FIG. 2(c). It can be 

20 observed from FIGs. 2(d) and 2(f), that while the interior (homogeneous) regions were 
(visually) equally smooth in the two images, contrast at the boundaries is significantly 
improved in FIG. 2(f). A similar phenomenon, although to a lesser degree, can be observed 
between FIGs. 2(c) and (e). 

It is xmreasonable to expect that a filtering method will retrieve structures that are lost 

25 by any degradation in the imaging process. Therefore, some finer details seen in FIG. 2(a) 
were lost in FIG. 2(f). Given the degraded phantom of FIG. 2(b) as the starting scene, it 
appeared in the filtered scene in FIG. 2(f) that someone had carefully smoothed the scene while 
preserving the structures and boundaries as much as possible. However, any such action would 
require at least a rough knowledge of the object size at various points, which is exactly what is 

30 provided by scale. 
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Phantom Example 1.2. 

The second phantom example was also a 2D phantom containing different geometric 
objects with varying shape, size and orientation. The entire phantom scene domain was divided 
into six square regions of equal size, three of which are at the top, and the other three are at the 
5 bottom. The initial phantom consisted of only two different intensity regions representing the 
objects (brighter) and the background (darker). This was generated by filling each of the six 
square regions with geometric objects of different size, shape and orientation. 

From this initial binary scene, 25 phantom scenes were generated by blurring (using a 
2D Gaussian kemel) the scene at 5 different degrees, and by adding a zero-mean Gaussian- 

10 correlated noise component at 5 levels. Each of the 25 scenes was further modified by adding a 
slow varying (ramp) background component which changed from 0 at the first column to 200 at 
the last column. Each of the resulting 25 scenes was filtered using the three filtering methods. 
One of the final phantom scenes, corresponding to a low level of blxirring and a medium level 
of noise, is shown in FIG. 6(a). The size of the phantom scene was 800x600. 

1 5 Due to the large size of the phantom, subtle details are not visible in the display of FIG. 

6(a). The phantom scene, as well as the results using the three filtering methods, are shown in 
FIG. 7 over each of the six ROIs marked in FIG. 6(a). FIG. 7(a) depicts scenes within the six 
ROIs taken from the unfiltered phantom scene, while FIGs. 7(b)-7(d) depict the corresponding 
filtered scenes that resulted from anisotropic diffusion ((shown in FIG. 7(b)), from scale-based 

20 neighborhood averaging (shown in FIG. 7(c)), and from scale-based diffusion (shown in FIG. 
7(d)). 

FIG. 6 also compares the three filtering methods via the intensity profile taken along a 
horizontal line, as shown in FIG. 6(a). FIG. 6(b) shows in magnification the profile over a 
small region indicated at the bottom of FIG. 6(a) (see rectangle outlined in black), and FIGs. 

25 6(c)-(f) show four profiles corresponding to the same region in four different scenes - the initial 
phantom scene after adding only the ramp component (but not adding noise and blurring), and 
the filtered scenes obtained using anisotropic diffusion (shown in FIG. 6(d)), using scale-based 
neighborhood averaging (shown in FIG. 6(e)), and using scale-based diffusion (shown in FIG. 
6(f)), respectively. In FIG. 6(d) (anisotropic filtering), the identities of the five (5) narrow, 

30 vertical columns, each representing a small square (see FIG. 7, 4th column), were lost; whereas 
by comparison, they were preserved by the scale-based methods (FIGs. 6(e) and 6(f)). 
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Meanwhile, filtering over each piecewise homogeneous region was also improved using 
scale-based diffusion, as compared with anisotropic diffusion. This was also visible in the 
background region shown in the ROI scenes. The first and the sixth columns in FIG. 7 were 
used to demonstrate the behavior of the three methods in homogeneous regions, while other 
5 columns demonstrate their behaviors in regions containing fine structures. Although all 
methods preserve large structures, the scale-based methods produced significantly improved 
results for finer structures. 

Example 2 - Quantitative Comparisons. 

10 The quantitative comparison utilized the scenes depicted in FIGs. 6 and 7. Let Ct = (C, 

fb) denote the original binary phantom scene (before adding artifacts), and IctE^ {P\ Q={Q 
fi), l<i<25} denote the set of 25 phantom scenes generated from as described in Example 
1, Phantom Example 1.2. Let let E^i = {C,/ 1 C,/ = (CJxiX l<i< 25} denote the set of scenes 
produced by applying the filtering method x e {ad, sa, sd} to the scenes in E, where ad, sa, sd 

15 correspond to anisotropic diffusion, scale-based neighborhood averaging and scale-based 
diffusion, respectively. Over the scene domain C, 100 ROIs, (f\ 1 < r < 100 were randomly 
selected, each of size 18x18, from homogeneous regions, such that all the pixels in a ROI had 
the same intensity (object or background) in the initial phantom scene Of the 100 ROIs, 60 
were firom background regions, and 40 were from object regions. Additionally, 50 ROIs, Cf\ 1 

20 < / < 50, were selected firom regions containing boundaries. Each Cf^ was selected such that it 
contained pixels firom both the object and backgroimd regions. Let 0/ and 5/ denote the set of 
pixels in (f^ from the object and background regions, respectively. 

For X e {ad, sa, sd} 1 < z < 25, and 1 <y < 100, "residual noise" was defined as 
the standard deviation of pixel intensities fxi(c) over each homogeneous ROI C^^ in a filtered 

25 scene resulting from method x for the phantom scene at the level of blurring and noise given by 
the index /. 

Tables 1, 2, 3 Ust in each cell, the mean and the standard deviation of RNl^^ , iW^^- , and 

RNI^. , respectively, at each degree of blurring and noise for the 100 ROIs. The degrees of 
blurring and noise increase along columns and rows, respectively, from the upper left comer. 
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Table L The mean (first entry) and standard deviation (second entry) of RN^^^^ values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
from left to right while blurring increases along columns from top to bottom. 





Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise 5 


tSlur 1 


o.Uio 


Q /lie 


iu.y4o 


1 "? A^l 
Lj.'tJ 1 






0.966 


1.278 


1.632 


1.984 


2.341 


Blur 2 


6.076 


8.498 


10.995 


13.506 


16.056 




0.985 


1.288 


1.637 


1.983 


2.338 


Blur 3 


6.035 


8.503 


11.086 


13.683 


16.606 




0.973 


1.298 


1.658 


2.015 


2.395 


Blur 4 


6.008 


8.512 


11.115 


13.691 


16.618 




0.965 


1.300 


1.657 


2.018 


2.401 


Blur 5 


5.987 


8.502 


11.111 


13.694 


16.625 




0.963 


1.299 


1.658 


2.021 


2.401 



Table 2 . The mean (first entry) and standard deviation (second entry) of RN^^. values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
from left to right, while blurring increases along columns from top to bottom. 





Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise 5 


Blur I 


4.655 


6.493 


8.813 


11.118 


13.490 




1.350 


1.491 


1.797 


2.157 


2.562 


Blur 2 


4.683 


6.661 


8.938 


11.244 


13.602 




1.339 


1.611 


1.893 


2.226 


2.598 


Blur 3 


4.603 


6.684 


9.282 


11.896 


15.586 




1.214 


1.638 


1.976 


2.354 


2.728 
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Blur 4 


4.580 


6.806 


9.415 


11.913 


15.604 




1.170 


1.543 


1.932 


2.371 


2.743 


Blur 5 


4.578 


6.795 


9.404 


11.920 


15.613 




1.161 


1.521 


1.919 


2.376 


2.747 



Table 3 . The mean (first entry) and standard deviation (second entry) of iW^^. values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
from left to right, while blurring increases along columns from top to bottom. 





Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise 5 


Blur I 


4.204 


5.920 


7.954 


10.072 


12.334 




0.923 


1.250 


1.631 


2.038 


2.530 


Blur! 


4.223 


5.966 


7.979 


10.114 


12.401 




0.924 


1.269 


1.646 


2.051 


2.529 


Blur 3 


4.238 


5.987 


8.254 


10.778 


14.774 




0.924 


1.276 


1.715 


2.235 


2.855 


Blur 4 


4.239 


6.100 


8.417 


10.785 


14.785 




0.924 


1.298 


1.747 


2.241 


2.854 


Blur 5 


4.240 


6.102 


8.420 


10.787 


14.797 




0.924 


1.298 


1.748 


2.240 


2.857 



From Tables 1, 2, and 3, it can be seen that at every level of blurring and noise, mean 
residual noise is maximal for anisotropic diffusion, and it is minimal for scale-based diffusion. 
Moreover, at each level of blurring and noise, three paired r-tests (Rosner, Fundamentals of 
Biostatistics , New York, NY, Duxbury Press, 1995, each over 100 pairs (iW^,- RN'y-),x^y,x, 
y e {ad, sa, sd}, and 1 < 100, of 7?/^ data produced by two different filtering methods, were 
conducted under the null hypothesis that there is no statistically significant difference between 
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the two methods. The hypothesis was rejected at a very high level of confidence (p < 10'^) in 
all tests between a scale-based method and anisotropic diffusion. This confidence level was 
less than 0.3x10"^ for the r-test between the two scale-based methods. 

Obviously, a measure of residual noise alone is not sufficient to characterize the 
effectiveness of a filtering method. For this reason, a measure of "relative contrast" was 
defined, RCi^ forx e {ad. sa, sd},l<i< 25, 1 <j < 50 as follows: 



RCi = 



(Equation 18) 



10 wherein M^' and cr^ denote the mean and standard deviation, respectively, of pixel intensities 
f^i(c) over the object region Oj and and crfj denote similar entities over the background 
region Bj. Notably, RC gives a measure of object-to-background contrast relative to residual 
noise in each region. 

Tables 4, 5, and 6 Ust in each cell the mean and the standard deviation of i?C , RC , 
1 5 and i?C^^,. , respectively, at each degree of blurring and noise for the 50 ROIs. Notably, a 
higher value of RC signifies more accurate filtering. 

Table 4. The mean (first entry) and standard deviation (second entry) of RC values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
fi-om left to right while blurring increases along columns fi:om top to bottom. 

20 





Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise 5 


Blurl 


7.983 


6.448 


5.185 


4.048 


3.240 




3.400 


2.076 


1.317 


0.772 


0.478 


Blur 2 


4.764 


3.691 


2.933 


2.404 


2.059 




2.048 


0.908 


0.307 


0.345 


0.491 


Blur 3 


2.286 


2.002 


1.784 


1.674 


1.609 




0.348 


0.561 


0.651 


0.682 


0.641 
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Blur 4 


1.833 


1.735 


1.611 


1.475 


1.418 




0.673 


0.696 


0.726 


0.751 


0.696 


Blur 5 


1.638 


1.561 


1.468 


1.365 


1.307 




0.755 


0.749 


0.753 


0.758 


0.700 



Table 5. The mean (first entry) and standard deviation (second entry) of RC^^^ values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
fi-om left to right while blurring increases along columns from top to bottom. 





Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise 5 


Blur I 


7.755 


6.340 


5.372 


4.627 


4.010 




2.645 


1.455 


0.860 


0.556 


0.412 


Blur 2 


5.176 


4.262 


3.635 


3.140 


2.705 




1.942 


1.138 


0.692 


0.429 


0.230 


Blur 3 


3.311 


2.679 


2.184 


2.002 


1.860 




0.919 


0.289 


0.244 


0.344 


0.369 


Blur 4 


2.390 


2.116 


1.878 


1.675 


1.590 




0.187 


0.266 


0.402 


0.485 


0.462 


Blur 5 


1.881 


1.780 


1.638 


1.501 


1.426 




0.396 


0.434 


0.495 


0.528 


0.502 



Table 6. The mean (first entry) and standard deviation (second entry) of RC^^- values 
are shown in each cell for different blurring and noise conditions. Noise increases along rows 
from left to right while blurring increases along columns from top to bottom. 
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Noise 1 


Noise 2 


Noise 3 


Noise 4 


Noise 5 


Blur 1 


8.033 


6.569 


5.563 


4.808 


4.181 




3.044 


1.721 


1.019 


0.662 


0.485 


Blur 2 


5.341 


4.386 


3.738 


3.250 


2.847 




2.238 


1.328 


0.809 


0.503 


0.294 


Blur 3 


3.503 


2.923 


2.321 


2.084 


1.913 




1.075 


0.486 


0.246 


0.381 


0.406 


Blur 4 


2.690 


2.226 


1.905 


1.673 


1.583 




0.423 


0.237 


0.445 


0.585 


0.530 




1.902 


1.766 


1.612 


1.481 


1.412 


Blur 5 


0.418 


0.492 


0.577 


0.627 


0.566 



From Tables 4 and 6 it can be seen that, at every level of blurring and noise, the mean 
relative contrast for the scale-based filtering method is higher than that for anisotropic 

5 diffusion. At each level of blurring and noise, three paired /-tests each of 50 pairs of (i?C^. , 
RC^y. \x^y,x,y e {ad, sa, sd), and 1 <7 < 50, of RC data produced by two different filtering 
methods were conducted under the null hypothesis that there is no statistically significant 
difference between the two methods. 

For the Mest between anisotropic diffusion and scale-based diffusion, the hypothesis 

10 was rejected at a very high level of confidence {p < O.SxlO'"^) at every level of blurring and 
noise, except for two (2) cases indicated by {Blur 1, Noise 1) and {Blur 1, Noise 2\ in which 
two cases the results were not statistically significant. For the Mest between anisotropic 
diffiision and scale-based neighborhood averaging, the hypothesis was rejected at a very high 
level of confidence (p < 0.28x10"^) at every level of blurring and noise, except for three (3) 

1 5 cases, which were at minimum blurring and the three lowest noise levels, and in which three 
cases the results were not statistically significant. For the ^test between the two scale-based 
methods, scale-based neighborhood averaging and scale-based diffiision, the hypothesis was 
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rejected at a very high level of confidence (p < 0.6x1 0'"^) at every level of blurring and noise, 
except for six (6) cases, which were at maximum blurring and the four highest noise levels and 
at second maximum blurring and the two highest noise levels, and in which six cases the results 
were not statistically significant. 
5 In conclusion, qualitative comparisons based on geometric phantoms generated under a 

range of conditions of blurring, noise and background variation, as well as patient MR images, 
indicate significant improvements and clear superiority using the scale-based filtering methods 
of the invention over the prior art anisotropic diffusion methods, for preserving boundary 
sharpness and fine details while suppressing noise. As between the two scale-based methods, 

10 results using scale-based diffiision have been found to be superior, except at extremely low 
blurring and noise where the differences were not statistically significant. Both scale-based 
methods embodying the invention use region-constancy based on local homogeneity to effect 
filtering. However, scale-based neighborhood averaging uses it in a hard fashion, while scale- 
based diffiasion uses it in a fiizzy manner. Because of this difference in the basic strategy, 

15 results using scale-based diffiision were considered to be more visually pleasing than those 
using scale-based averaging. Since the same trend was observed in the quantitative 
experiments, scale-based diffiision was preferred over scale-based neighborhood averaging, 
although the latter is computationally the least expensive among the three methods. 

Each and every patent, patent application and publication that is cited in the foregoing 

20 specification is herein incorporated by reference in its entirety. 

While the foregoing specification has been described with regard to certain preferred 
embodiments, and many details have been set forth for the purpose of illustration, it will be 
apparent to those skilled in the art that the invention may be subject to various modifications 
and additional embodiments, and that certain of the details described herein can be varied 

25 considerably without departing from the spirit and scope of the invention. Such modifications, 
equivalent variations and additional embodiments are also intended to fall within the scope of 
the appended claims. 



389797-5 



31 



M-2242 



